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Summary. A path-integral approach to lattice polarons is developed. The method 
is based on exact analytical elimination of phonons and subsequent Monte Carlo sim- 
ulation of self-interacting fermions. The analytical basis of the method is presented 
with emphasis on visualization of polaron effects, which path integrals provide. Nu- 
merical results on the polaron energy, mass, spectrum and density of states are given 
for short-range and long-range electron-phonon interactions. It is shown that cer- 
tain long-range interactions significantly reduce the polaron mass, and anisotropic 
interactions enhance polaron anisotropy. The isotope effect on the polaron mass and 
spectrum is discussed. A path-integral approach to the Jahn- Teller polaron is devel- 
oped. Extensions of the method to lattice bipolarons and to more complex polaron 
models are outlined. 



1 Introduction 

The polaron problem was one of the first applications of path integrals (PI). 
Just a few years after the introduction of the new technique into quantum 
mechanics [H [2] and quantum statistics [3] Feynman published his seminal 
paper [3] on the theory of Frohlich polaron [S], which set a stage for the 
development of polaron physics for the next half a century. The main feature 
of the statistical PI is an extra dimension that transforms each point-like 
particle into a one-dimensional line, or path. The extra dimension, to be called 
here the imaginary time r, extends between and /3 = (fcsT)^^ where T is 
the absolute temperature at equilibrium. As a result, a quantum-mechanical 
system is mapped onto a purely classical system in one extra dimension. This 
enables visualization of quantum-mechanical objects, which is an instructive 
and useful feature. The quantum mechanical properties are fully recovered by 
considering an ensemble of paths, each contributing its own statistical weight 
into the system's partition function. 

Having formulated the polaron problem in PI terms, Feynman made an- 
other key advance. He showed that if the ionic coordinates entered the model 
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only quadratically (a free phonon field) and linearly (electron-phonon inter- 
action in linear approximation) then the infinite-dimensional integration over 
ionic configurations could be performed analytically and exactly. As a result, 
all bosonic degrees of freedom (an infinite number of them) are eliminated 
in favour of just one fermionic degree of freedom. Phonons remain in the 
system as a retarded self-interaction of the electron. In the PI language, the 
statistical weight of each path is exponential in its Euclidian action and the 
latter contains a double integral over imaginary time as opposed to a single 
integral in cases of ordinary instantaneous interactions. Different segments of 
the fcrmion path "feel" each other if they are separated in imaginary time by 
less than the inverse phonon frequency. As a result, the electron path stiffens 
which leads to an enhanced effective mass and other characteristic effects as 
detailed later in the chapter. 

The resulting electron path integral could not be calculated analytically 
due to the complex nature of the retarded self-interaction. Feynman resolved 
the difficulty by employing an original variational principle, in which the ex- 
act polaron action was replaced by an approximate, but exactly solvable, 
quadratic action. That led to a remarkably accurate top bound for the energy 
of the Frohlich polaron (see [6l[7l[8]), which was only marginally improved by 
subsequent generalizations of the method [HI UHl [HI [13 [H] ■ The PI method 
became a workhorse of the Frohlich polaron research for decades and inspired 
numerous extensions which included polaron mobility |141 1151 116j and optical 
conductivity [l7] , polaron in a magnetic field [HI [19] , large bipolaron [20] , 
and others. Feynman's method also inspired a Fourier Monte Carlo method, 
in which only the difference between the exact and variational energy was 
calculated numerically [211 [Hj. Recently, the variational PI treatment was 
extended to a many-polaron system [23| . 

In the parallel development of the small-polaron physics, Feynman's re- 
markable reduction had not been utilized for almost thirty years. Although 
the phonon integration could be performed as well the self-interacting elec- 
tron PI could not be calculated. The problem was that on the lattice even a 
free particle possessed a non-quadratic action. Because such a path integral 
could not be done analytically, it could not serve as a trial variational action. 
The situation changed in 1982 when De Raedt and Lagendijk (DRL) observed 
[24] that the electron PI could instead be evaluated numerically. Metropolis 
Monte Carlo [55] was ideally suited for the task because the polaron action 
was purely real, and as such resulted in a positive-definite statistical weight of 
the path. Using this approach, DRL obtained a number of nice results on the 
Holstein polaron: confirmed formation of a self-trapped state with increasing 
strength of the electron-phonon interaction, observed that the crossover gets 
sharper with decreasing phonon frequency and increasing lattice dimensional- 
ity [24] , and that the critical coupling goes down as a dispersive phonon mode 
softens [27]. These results were reviewed in [25]. DRL also provided the first 
Monte Carlo analysis of the Holstein bipolaron [28] . 
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The works of DRL were an important step forward in the apphcation of 
PI to lattice polarons. Their method could handle infinite crystal lattices of 
arbitrary symmetry and dimensionality, dispersive phonons, and long-range 
electron-phonon interactions. At the same time, it was still limited to ther- 
modynamical properties such as energy, specific heat, and static correlation 
functions. In addition, it suffered from a systematic error caused by finite 
discretization of imaginary time (the Trotter slicing). These limitations were 
removed in 1997-1998. Firstly, it was shown how open boundary conditions in 
imaginary time could enable direct calculation of the polaron effective mass 
[29] and even the entire spectrum and density of states [30] . The open bound- 
ary conditions allow projection of the partition function on states with definite 
polaron momenta. This is a particular manifestation of the general projection 
technique in the presence of a global symmetry. Systematic application of this 
principle makes it possible to separate states of different symmetries in many 
interesting situations, most notably bipolarons of different parity and orbital 
symmetry. More about the projection technique is in Section [2] Secondly, the 
polaron Monte Carlo method was formulated in continuous imaginary time 
[31j . which completely eliminated the Trotter slicing and tremendously im- 
proved the computational efficiency of the method. This will be described in 
detail in Section [3] As a result, the versatility of DRL's method was enhanced 
by better computational efficiency and by a number of new polaron and bipo- 
laron properties that could be computed with path integrals. The method was 
further developed in [32 [Ml IM 133 [Ml EZ] |3H1 [39] , the content of which will 
be described later in the Chapter. 

The path-integral Quantum Monte Carlo (PIQMC) with phonon integra- 
tion is only one from an impressive list of numerical methods developed for 
polaron models in the last three decades. There exist several other QMC 
techniques, in particular the path-integral QMC without phonon integration 
[40l[4l], Fourier QMC [211 [22], diagrammatic QMC [H [H] [42l [43l H [45] and 
Lang-Firsov QMC [13 [HI [IS]. Non-QMC classes of methods include exact di- 
agonahzation [13 [501 [HI [SI [Ml [Si], variational calculations [TOl [551 [551 [57 1 1551 
[59l l60l [611 [62l [63] and the density-matrix renormalization group [gH [65l [66] . 
Despite proliferation of methods, most of them have been applied to the two 
major polaron models: the ionic crystal model of Frohlich [5] and molecular 
crystal model of Holstein [67]. (The notable exceptions are the Jahn- Teller 
(bi)polaron [Ml [50] and Peierls instabilities [HI [BSl US].) Due to its versatil- 
ity, PIQMC is uniquely positioned to study many other models, for example 
long-range and anisotropic electron-phonon interactions. In addition, it pro- 
vides some (bi)polaron properties that are difficult to obtain by other methods 
such as the density of states. In this way, several physically interesting results 
have been obtained that include, in particular, the light polaron mass in the 
case of long-range electron-phonon interactions [32] , the enhancement of po- 
laron anisotropy by electron-phonon interaction [33] . formation of a peak in 
the polaron density of states [Ml [Ml [37] , the isotope effect on the polaron 
spectrum and density of states [35, and the "superlight" crab- like bipolaron 
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[701 [m 133 ■ These and other results obtained by the PIQMC method in the 
last decade will be reviewed in Section 21 Various extensions of PIQMC and 
its prospects are given in Section [51 

2 Projected partition functions 

It is commonly believed that the statistical path integral can provide informa- 
tion only on the ground state of a quantum mechanical system, and in general 
this is true. However, when the system possesses a global symmetry suitable 
projection operators can project the configurations on sectors of the Hilbert 
space corresponding to different irreducible representations of the symmetry 
group. That enables access to the lowest states within each sector and pro- 
vides valuable information about system's excitations. This strategy proves 
very useful in the path-integral studies of polaron models. It enables calcula- 
tion of the polaron mass, spectrum, density of states, bipolaron singlet-triplet 
splitting and other properties. The above considerations are formalized as 
follows. The full thermodynamic partition function is a trace of the density 
matrix: 

Z = ^(R|e-^^|R), (1) 

R 

where H is the Hamiltonian and |R) is a complete set of basis states in 
the real-space representation. If there is a global symmetry group G with a 
set of irreducible representations J7, it is meaningful to compose a projected 
partition function Gjj which is a trace over the {/-sector of the Hilbert space 
only: 

Zu = Y.{Ku\e-^"\B.u) = Y.{no\je-P"Ou\R) , (2) 

R(7 R 

where operator Ou generates a basis state |Rc/) from an arbitrary state |R). 
In the low-temperature limit /3 — > oo, the partition function is dominated by 
the lowest {/-eigenvalue, Zjj exp {—(3Eu). Thus Ejj can be found by taking 
the logarithm of Zjj in the low-temperature limit. This is particularly efficient 
if the first excited state in U is separated from Ejj by a finite gap. 

The most important application of this technique is the formula for the 
(bi)polaron effective mass. In this case, irreducible representations U are 
labeled by the total momentum K and the projection operator is Ok = 
^^giKrjT^^ where is the shift operator by a lattice vector r. The K- 
projected partition function is then ^30] 

Z^ = Y, e*K^'-(R + AY\e-f'"\R) = ^ e^^^^^Z^r ■ (3) 

Ar Ar 

Here R-|-Z\r denotes a many-particle configuration R, which is parallel-shifted 
as a whole by a lattice vector Ar. The partition function Z^r is a "shifted 
trace" of the density matrix: it connects R not with R but with R shifted by 
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Ar. It follows from the above equation that Zk and Z/^^ satisfy a Fourier-type 
relationship (for a more detailed derivation, of. Refs. [3D1[72])- The partition 
function Z^^ is formulated in real space and as such can be represented by 
a real-space path integral. The partition function Zk is diagonal in momen- 
tum space and therefore contains information about the variation of system's 
properties with K. The Fourier theorem ^ is key to the ability to infer the 
(bi)polaron spectrum and mass from the path integral. 

Dividing Zn by Zk=o one obtains in the zero temperature limit 

The ratio of two sums over Ar represents the mean value of exp (iKZ\r) over 
an ensemble of paths in which the two end-point many-body configurations 
are identical up to an arbitrary parallel shift. In the QMC language, simula- 
tions need to be performed with open boundary conditions in imaginary time. 
There is a certain parallel with the uncertainty principle: fixing the boundary 
conditions in imaginary time (by making them periodic) results in a mix of all 
possible momenta. Conversely, relaxing the boundary conditions by mixing in 
all possible shifts Ar allows extraction of information on the given K-sector of 
the Hilbert space. Resolving the last equation with respect to E-k, one obtains 
an estimator for the (bi)polaron spectrum 

Ek- Eo = - lim i ln(cos KZ\r)shift ■ (5) 

P^oc p 

The left-hand-side of this equation is a constant that depends on the pa- 
rameters of the model being simulated. Therefore, the average cosine on the 
opposite side must decrease exponentially with /? to compensate the growing 
denominator. At some f3 it becomes too small to be measured reliably with 
the available statistics of a Monte Carlo run. This is an incarnation of the 
infamous "sign problem" that plagues many QMC algorithms. In physical 
terms, at very low temperatures the probability to access an excited state is 
exponentially low due to the Boltzmann weight factor. As a result, the QMC 
process samples relevant configurations exponentially rarely. These considera- 
tions put the low limit on the temperature of the simulation (high limit on (3). 
At the same time, the condition that Zk is dominated by a single eigenstate 
in the K sector, puts a high limit on temperature. Thus there is an interval of 
temperatures at which a meaningful calculation of (bi)polaron spectrum can 
be performed. If the expected bandwidth is too large the temperature inter- 
val shrinks to zero, which renders the calculation impossible. In general, the 
(bi)polaron spectrum can be computed by this method if the total bandwidth 
is smaller than the phonon frequency. 

Expanding ^ for small momenta, one derives an estimator for the ii-th. 
component of the inverse effective mass 

^ = lim . (6) 
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Since in the present formulation the end points of the path are not tied to- 
gether path evolution between r = and t = [3 can be regarded as diffusion 
of the top end with respect to the bottom end over a diffusion time /3. Ac- 
cording to the Einstein relation, the mean square displacement of a diffusing 
particle is proportional to time. Then ([6]) implies that the inverse effective 
mass is proportional to the diffusion coefficient D of the path in imaginary 
time [29l [3TJ [72] • The combination h(3 has the dimensionality of time, so the 
inverse mass can be written as 

^ = ¥- (7) 



This formula works for both polarons and bipolarons [441 [39]. It can be viewed 
as a fluctuation-dissipation relation because it equates a dynamical character- 
istic, the mass, with an equilibrium thermodynamic property D. In addition, it 
provides a nice visualization of the main polaron effect: effective mass increase 
caused by an electron-phonon interaction. As will be shown in the next sec- 
tion, phonon integration results in correlations between distant parts of the 
imaginary-time polaron path. This increases the statistical weight of paths 
with straight segments. The average polaron path stiffens, which translates 
in a reduced diffusion coefficient D and, according to 0, enhanced effec- 
tive mass. Note that these considerations apply equally well to other types of 
composite particles whose mass originates from interaction. The most notable 
examples are defects in quantum liquids [73] and the hadrons of quantum 
chromodynamics . 

Another important application of the projection technique is the singlet- 
triplet splitting of the bipolaron, or, more generally, of a two-fermion bound 
state. Coming back to the projected partition function Zjj, the relevant 
projection operators are identified as Os = I + X for singlet states and 
Ot = I — X for triplet states. Here / is the identity operator and X is 
the fermion exchange operator, X|rir2{^}) = |r2ri{^}), where {^} are the 
ionic displacements. Upon substitution in ^ one finds 

Z^-^^ J2 (rir2U}|e-^^|rir2U}>±(r2ri{a|e-''^|rir2{a). (8) 

rir2{4} 

Monte Carlo simulation proceeds over an ensemble of paths that have iden- 
tical end configurations up to exchange of the top ends of the two fermion 
trajectories. For the singlet, both types of configurations contribute a phase 
factor (-1-1) whereas for the triplet the direct paths contribute (-1-1) and ex- 
changed paths contribute (—1). Forming the ratio of the triplet and singlet 
partition functions one obtains an estimator for the S — T energy split 

-E^ = - ]im (9) 

Here (—1)^ = ±1 for the direct and exchanged two-fermion paths. If the ends 
of the two paths coincide then (—1)^ = 0. All said above about the sign 
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problem and the role of temperature in calculating the average cosine applies 
to this expression as well. 

Equation ^ provides the energy difference between the lowest triplet and 
singlet states. It is possible to compute the effective mass and the entire spec- 
trum of the triplet bipolaron. To this end, one needs to combine the two 
symmetries considered above: the translational and exchange. This naturally 
leads to a path ensemble with the end configurations differing by any combi- 
nation of parallel shifts and exchanges. Without repeating the derivation, the 
final expressions are 

,T _ 1 (ZlZl\ _ 1 ,^ ((-l)^cos(KZ\r)),hift,c 



for the triplet spectrum and 

1 1 ((-l)^(Z\r02)shift,, 



13%^ ((-l)^)shift,, 



(11) 



for the triplet mass. For the singlet spectrum and mass, ([5|) and ([6]) are still 
valid, except that the boundary conditions must be changed to "shift and 
exchange" . The same applies to the singlet-triplet estimator ([9]) . 



3 Continuous-time path-integral quantum Monte Carlo 
method 

The main appeal of quantum Monte Carlo (QMC) methods is the ability to 
calculate physical properties without approximations. A quantum mechanical 
system is directly simulated taking into account all the details of particle dy- 
namics and inter-particle interaction. Physically interesting observables can 
be calculated without bias while the statistical errors in general decrease with 
increasing simulation time. It is said therefore that QMC provides "numeri- 
cally exact" values for the observables. A number of comprehensive reviews 
on QMC exist [H [IS US [76l [77l [TS]. Many problems of the early QMC 
methods, such as the finite-size effects, finite time-step effects, and critical 
slowing down, have been resolved with the development of novel algorithms 
and increasing computing power. The sign problem, that is the non-positive- 
definiteness of the statistical weight of basis configurations, remains the only 
fundamental problem. In real systems and models without a sign problem (for 
example liquid and solid He"* [HI ES HOl [H HI] ) , QMC works beautifully and 
provides with accurate and valuable information about the thermodynamics 
and sometimes real-time dynamics. 

One polaron is free from the sign-problem difficulties. Phonons are bosons 
and as such do not lead to a fermion sign problem. And as long as there is 
only one polaron in the system statistics does not matter. This is the main 
reason behind the success of the QMC approach to the polaron. The ground 
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state energy, the effective mass, isotope exponents on mass, the number of 
excited phonons, static correlation functions and other quantities can be cal- 
culated without any approximations. The bipolaron ground state can also be 
investigated without a sign problem as long as the bipolaron is a spin singlet 
with a symmetric spatial wave function. Thus the power of QMC can be (and 
has been) applied to bipolarons of various kinds and to pairs of distinguish- 
able particles such as the exciton. For many-polaron systems the fermion sign 
problem fully manifests itself, which limits the applicability of QMC meth- 
ods. More about this will be said in Section [S] In this Section the basics of 
the continuous-time path-integral quantum Monte Carlo (PIQMC) method 
for lattice polarons are described. 

The starting point is the shift partition function Z/ir, see ([3]), and the 
electron-phonon (c-ph) Hamiltonian 

i/^-t^ci,c„-^/.(n)ci,c„U + E("|^^ + ^e) . (12) 

nn' nm m ^ ™ ' 

The Hamiltonian is written in mixed representation, Cn being fermionic an- 
nihilation operators and ion displacements. Index n denotes the spatial 
location of an electron Wannier orbital whereas m denotes that of an ion dis- 
placement. In general, n 7^ m even if they belong to the same lattice unit cell. 
For simplicity, the electron kinetic energy (the first term of H) is taken in the 
nearest-neighbor-hopping approximation with amplitude — t, and the lattice 
energy (third term in H) in the independent Einstein oscillator approximation 
with mass M and frequency uj. Neither approximation is necessary. Section O 
explains how PIQMC deals with more complex forms of kinetic energy and 
phonon spectrum. The most interesting part is the c-ph interaction (second 
term in H) which is written in the "density-displacement" form. The quantity 
/m(n) is the force with which an electron n acts on the ion coordinate m. 
Asymmetric notation emphasizes that the force /m is an attribute of a given 
oscillator, while the argument n is a dynamic variable indicating the current 
position of the electrons interacting with m. No constraints are imposed on 
the functional form of /, thus allowing studies of long-range e-ph interactions. 
The commonly used Holstein model [67, corresponds to a localized force func- 
tion /in(n) = KSnm- 

3.1 Handling kinetic energy in continuous time 

Handling the kinetic energy on a lattice possesses a challenge for Monte Carlo. 
To see this consider just the kinetic part of H and one free particle. By 
introducing multiple resolutions of identity, / = J2r \^){^\' the shift partition 
function is developed into a multidimensional sum 

ZAr= (ri+/ir|e-^^^-|rM)...(r3|e-^^^-|r2)(r2|e-^^^-|ri), 

rir2...rA/ 

(13) 
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where At — [3 / L is the time step and L + 1 is the number of time slices. Each 
matrix element from the product is expanded to the first power in At: 

(r,+ile-^^^-|r,) = 5r,,„r, + Mr ^ <5r,.^„ + 0{At^) , (14) 

1 

where 1 runs over the nearest neighbors. Expanding the product, the partition 
function becomes a sum of a large number of terms, each of which represents 
a particular path of the electron in imaginary time r(r). The paths consist 
of two different building blocks: straight segments and kinks, which originate 
from the first and second terms in p^ . respectively. On straight segments 
the electron position does not change, whereas each kink changes the electron 
coordinates by vector 1 depending on the kink type. The statistical weight of 
a path with Nk kinks is 

Wn, = (Mr) (Mr) . . . {tAr) = (Mt)^'^ . (15) 

Nk times 

Now consider two paths, D and D', that differ by D' having one extra kink 
of type 1 at imaginary time r. At small At, D' has a vanishingly small sta- 
tistical weight compared to D. As a result, no meaningful detailed balance 
between individual paths with different number of kinks seems to be possible 
in the continuous-time limit. The resolution of this difficulty comes from the 
observation that there are infinitely many paths with A^^ + 1 kinks than with 
N}^ kinks. The small weight of individual path is compensated by the large 
number of paths. The combined weight of all paths with A'fc -I- 1 kinks is of 
the same order than the combined weight of all paths with A^j. kinks. Instead 
of detailed balance between individual paths one should seek detailed balance 
between path groups of similar combined weight. The paths D and D' are 
regarded as representatives of such groups. 

In accordance with the general ideas of the Metropolis Monte Carlo [26] . 
one has to compose an equation that equates the transition rate from D to 
D' to the reciprocal transition rate from D' to D. Both rates are products of 
three factors: the probability to have the original path, which is proportional 
to the path's weight W, the probability Q to propose the necessary change to 
the path, and the probability P to accept the change. The detailed balance 
reads: 

W^(L>)Q(add 1 at t)P{D D') = T^(7:>')Q(remove 1 at t)P{D' D) . (16) 

The probabilities to propose the addition and removal of the kink turn out 
to be two quite different quantities. In the direction D' D, the kink to 
be removed is chosen from a finite number of all other kinks. Therefore the 
probability that it is selected for removal is a finite number Qremovo(l; ''■)• In 
contrast, in going from D to D' the kink does not yet exist and the probability 
to propose its creation exactly at time r is proportional to the time interval 
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At. The probability to propose the kink's addition is an infinitesimal quantity: 
Q(add 1 at r) = qi^^d{^;T){AT), where (/add is the corresponding probability 
density. Substitution of these two results in the balance equation yields 

{tATf''q,M{l;T){Ar)P{D ^ D') = (MT)^''+ig„move(l; r)Pp' ^ D) . 

(17) 

The extra (At) from the addition probability on the left side exactly balances 
the smaller weight of the right side. The time step cancels out of the equa- 
tion leaving only finite factors which are well defined in the continuous-time 
limit. According to the standard Metropolis recipe the solution is given by 
the following acceptance rules 

P{D ^ JO = min (l ; ^' ^) \ , (ig) 

L '7add(l;T) J 

P(i?'^I?)=min(l; -^^i^fej-rj . (19) 

These expressions still leave a lot of freedom in specifying functions gadd and 
Qromovc- As an example, consider the simplest choice when both functions 
are constant. For the addition process this means that the time for the new 
kink of type 1, on top of the existing is chosen with equal probability 
within the [0, j3\ time interval, hence (;add(l; t) ~ For the removal process, 
the candidate kink for elimination is chosen from the existing Ni with equal 
probability, that is Qremove(l; ■'') = -^i"^- One must also take into account the 
fact that kinks can only be added to the path if A^i = 0, while added or 
removed if TVi > 1. The full set of update rules could be as follows, (i) Select 
the kink type 1, that is the hopping direction of the particle. The probability 
with which different 1 are selected may not be equal but must be non-zero, (ii) 
If the current path already has kinks of type 1, whether to add or remove a kink 
is proposed with, for example, equal probability of ^. If the path possesses no 
kinks, addition is the only option, and it must be chosen with probability 1. 
(iii) If addition is selected, the time for the new kink is chosen with probability 
density i. The update is accepted with probability min{l; {t(3)/{N\ + 1)} for 
TVi > 1 and with probability min{l; (t/3)/2} for A'^i = 1. (iv) If removal is 
selected, the candidate kink is chosen with equal probability from the existing 
Ni. The update is accepted with probability min{l; A^i/(</3)} for A^i > 2 and 
with probability min{l; 2/(i/3)} for iVi = 1. 

The described method can be readily generalized to next-nearest neighbor 
hopping and anisotropic hopping. It can also be generalized to multi-kink up- 
dates, as needed for example in simulating the bipolaron. To close this section 
one should note that a path can be interpreted as a time-space diagram, where 
kinks play the role of vertices and straight segments the role of propagators. 
Integration over the ensemble of fluctuating paths can then be understood 
within the general concept of Diagrammatic Monte Carlo [7l [83] . 
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3.2 Integration over phonons 

Analytic integration over phonons is a key ingredient of the polaron PIQMC 
method. Although an electron-phonon (e-ph) system can be simulated with- 
out phonon integration [40l [41] , it is the integration that makes the method 
powerful. Indeed, the integration reduces A^-|- 1 degrees of freedom to just one 
degree of freedom, which increases accuracy and reduces simulation time. The 
integration effectively eliminates all the distant parts of the system by folding 
the infinite number of ion coordinates into one retarded self-interaction func- 
tion. The simulation then proceeds on an infinite lattice with no finite-size 
effects. In addition, the integration makes it possible to simulate long-range 
e-ph interactions with the same efficiency as local interactions, which extends 
QMC capabilities towards more realistic models. 

Technical details of the phonon integration are somewhat cumbersome but 
well-documented in the literature. Here we will describe the main steps of the 
derivation and point out associated subtleties. 

Upon insertion of H into ([31), a path- integral representation for Z/\r is 
developed by introducing L time intervals and L — 1 intermediate sets of basis 
states. In the continuous time limit, L — > oo, the kinetic energy of the electron 
is treated as described in the previous section. Integration over the electron 
coordinates is converted to stochastic summation over an ensemble of paths 
which are sampled by adding and removing kinks. For each electron path, 
there is an internal path integration over the ionic coordinates characterized 
by the electron-phonon action ^c-ph 

ZAr= Vv{t) P^^(T)e^--'' Ml, (20) 

"'(r.O) "'({«„},0) 



2h^ 



d^ + E //m[r(r)]C„.(r)dT. 



(21) 

The path integral over ^m(''') is Gaussian and therefore can be calculated ana- 
lytically. Actually, the integration is performed in two steps. Firstly, path inte- 
gration is done with fixed but arbitrary boundary conditions at the two ends of 
the path. Secondly, a certain relationship between ^m(O) and Cm(/3) is assumed 
followed by an additional one-dimensional ordinary (but still Gaussian) inte- 
gration. In most PI studies, periodic boundary conditions ^m(/5) = Cm(0) are 
assumed. This leads to the classic Feynman formula for the oscillator action in 
the presence of an arbitrary time-dependent external force [84 . However, we 
have seen in the previous section that mass and spectrum calculation requires 
a partition function with the end configurations parallel-shifted with respect 
to each other. The shift of the ionic configuration must be the same as that of 
the fermion configuration, which implies ^m(/^) = 'fm-zir(O), as indicated in 
the last equation. The correlation between the electron and phonon boundary 
conditions is illustrated in Fig. [IJa). 
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Fig. 1. (a) Integration over ionic paths under shifted boundary conditions. Ionic 
displacements are shown as vertical bars at r = and t — 13. Notice how the 
pattern of displacements is correlated with the shift of the electron path. The end 
displacements of individual oscillators are not equal, fm(0) 7^ ^m(/3). (b) Calculation 
of the polaron action as a double integral over imaginary time. The kink times 
break the (tt') plane into a finite number of rectangles. Within each rectangle 
the electron coordinates are constant and the double integral can be calculated 
analytically for arbitrary [36]. After that the double integral reduces to a sum 
over the rectangles. Since the number of rectangles decreases at strong coupling, the 
algorithm gets faster at strong coupling. 



Due to an additional shift, the second-phase integration results in a cor- 
rection to the conventional "periodic" action. The full polaron action, upon 
summation over all oscillators, reads [3T1 [55] 

ApoMr)] = Ao + AA , (22) 

^o = Et^ /rd.d.- -^^-(f:';--'') /^[r(.)]/^[r(/)], (23) 
^4wM7o7o smh?iwf 

^ E - Cr^) , (24) 

m 

Bni= dr . , ^ „ /mir(T )J , (25) 

Jo smn nuip 

C^^f^r'^^^f^Hr')]. (26) 
Jo smn TiLop 

The polaron action is a functional of the electron path and contains all the in- 
formation about the ionic subsystem. The denominator of ([4]) finally becomes 
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Fig. 2. (a) A typical fluctuating path of a free particle on a one- dimensional chain 
(_Eo = — 2.00t, m* /rriQ = 1). The path has drifted to the left by four lattice sites, 
which contributes to the effective mass according to (|6]). (b) A typical path of a 
one-dimensional Holstein polaron at Tiuj/t — 2.0 and A = 2.0 (£"0 = —4.66(2) 
m* /mo = 3.06(3). The path has drifted to the right by just one site reflecting a 
larger mass. The thick bar is a kink update proposed by Monte Carlo. If accepted, 
the top part of the path is changed to the dashed line (top shift). Alternatively, the 
bottom part of the path can be shifted by one site to the left {bottom shift). 



^K=o = E^^r = ^ph E / ••• / (dr)^^t"'^e^-'[-Ml. (27) 

At Nk=0,l.,...° ° 

Here = [2 smh{^'h/3uj)]^^ (N is the number of lattice sites) is the partition 
function of free phonons. It is a multiplicative constant that cancels out from 
all statistical averages. The factor e"^p°' adds to the weight of each path on 
top of the "kinetic" contribution (tZ\r)^'=. As a result, the acceptance rules 
(fT8|) and (fT9|) are modified by the ratio of the respective factors for the paths 
D and D' 

P{D ^ D') = min (1 ; ^-Q™c(1;t) ^a,^,(d')-a,^,(d)\ ^ ^38) 

L gadd(l;T) J 

P[D' -^D) = min (1 ; ^"'^'^''^^ ^A,^.(D)-A,^,iD-)\ (39) 

t • C/rcmovc(l5 ^) J 

These acceptance rules together with the formula for the polaron action 
and path update rules described in the previous section fully define the Monte 
Carlo update process. It is illustrated in Fig.O Numerical efficiency of the al- 
gorithm critically depends on how easily the polaron action can be computed 
at every update. The task is greatly aided by the fact that the electron path 
is defined on a discrete lattice. In other words, it consists of a finite number of 
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straight segments. This has two important consequences. First, because within 
each segment the electron coordinate is independent of imaginary time r, the 
forces /m(r) can be taken outside the time integration. The rest of the inte- 
grand is an exphcit function of r and r'. This time integral can be calculated 
analytically for arbitrary end times of the segment. The double integral over 
time becomes a double sum over the segments, as illustrated in Fig.[T]Jb). The 
summand is a simple function of kink times; the explicit expressions are given 
in [36] . The smaller the number of segments, that is the straighter the path, 
the less time is required to compute ^poi. The electron-phonon interaction 
increases the statistical weight of straight paths and reduces the mean num- 
ber of segments. Thus, in contrast to most other numerical methods, PIQMC 
becomes faster at large electron-phonon couplings. 

Secondly, the (zj)-th term of the sum over segments contains the electron- 
ion forces in the combination 

<?(r,-r,) = ^/™(r,)/m(r,), (30) 

m 

where is the electron coordinate on i-th segment. The coefficients <P are 
defined on a discrete set of points and can be pre-computed for a sufficiently 
large range of coordinate separation. After that the simulation takes essen- 
tially the same amount of time for any form of electron-phonon interaction. 
This is why the PIQMC polaron method is as efhcient for long-range e-ph 
interactions as for the short-range Holstein model. 

We close this section by noting that the requirement of phonon integration 
with shifted boundary conditions equally applies to the continuous case. This 
does not seem to be the case with the original Feynman's calculation of the 
mass of the Frohlich polaron and its subsequent generalizations. This issue 
was thoroughly investigated in |85j . It was shown that the periodic boundary 
conditions in the phonon integration, while absolutely correct for the energy 
calculation, lead to infinite terms in the polaron action in the mass calcula- 
tion. Feynman's mass formula is obtained only if the terms are regarded as 
unphysical and thrown out by hand. In contrast, the shifted boundary condi- 
tions result in a self-consistent action and Feynman result is obtained without 
any complications. 

3.3 Observables 

The preceding sections explained how to organize an efficient Monte Carlo 
sampling process that simulates the electron kinetic energy while fully taking 
into account the ion dynamics and electron-phonon interaction. Note that no 
approximation has been made so far except for the quadratic form of the 
elastic energy of the crystal. In this section it will be explained how to extract 
useful physical information from the one-to-one correspondence between the 
quantum mechanical-polaron and an ensemble of self-interacting paths. A 
number of important estimators have already been derived in Section [2l In 
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particular, the effective mass is obtained from the mean-square end-to-end 
distribution, see ©■ The polaron spectrum relative to the ground state is given 
by (0) ■ It is essential that statistics for different En is collected in parallel so 
that the entire polaron spectrum is obtained in a single run. One consequence 
of this remarkable property is PIQMC's ability to compute the polaron density 
of states by simply histogramming _Ek at the end. At present, PIQMC is the 
only method capable of efficiently calculating the polaron density of states. 

The next important characteristic is the absolute polaron energy Eq. 
One way to obtain Eq is from the low-temperature limit of the internal 
energy U. The latter can be computed from the thermodynamic relation 
U — -~dh\ Z / dp. Since we are interested in the ground state, the correspond- 
ing partition function must be Zk=o = Hav^Av = Eshiftod paths 
Returning momentarily to discrete time, (3 = LAt, 

rj ^ l Eshift^ ^ _ 1 EshiftiiMZ = _^ 

^ EshiftW^ L EshiftW^ M\wdAT/^^,,- > 

The weight is a product of the kinetic contribution (iZir)^* and the phonon 
term e"^p°'. Substituting in the above expression, one obtains a ground-state 
energy estimator 

E.-^-(^ + '4^) ■ (32) 



L\At dAr/^^,,, \ (3 8(3 , 

The two terms here represent the kinetic and potential energy of the polaron, 
respectively. The kinetic energy is simply the mean number of kinks on the 
path. One can regard the kinks as "quanta" of the kinetic energy, each con- 
tributing the same amount ~(3^^. With increasing coupling the paths become 
"stiffer" , that is the mean number of kinks goes down. As a result, the polaron 
kinetic energy decreases by absolute value. Measuring the potential energy is 
harder for it is derived from the polaron action. Similarly to the latter, the 
potential energy estimator can be reduced to a double sum over the path's 
segments. Explicit expressions for the double sum terms are given in [36| . 

The number of excited phonons A'ph in the polaron cloud is another inter- 
esting characteristic. It measures the amount of energy stored in the lattice 
deformation around the electron. The polaron mass scales exponentially with 
A'ph so in addition the latter gives a rough estimate of the polaron mass. An 
estimator for A^ph can be derived by noting that the last term in the Hamil- 
tonian (fT^ can be rewritten as hu!{Nph + ^), where iVph is the operator of 
the total number of phonons. The mean value of A^ph is obtained by differen- 
tiating the free energy with respect to huj. To make sure that no contribution 

is received from the interaction term, the derivative must be taken with the 

1 



combination p/{Mu)) kept constant. Using F — — -ilnZ, the estimator is 



obtained as follows 



~ i3 d{huj) /2 - ^ \ dihuj) /2 / • ^^^^ 

^ ' — \ \ ' —/ shift 
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The subscript "shift" implies that the number of phonons corresponds to the 
ground state of the polaron (K = 0). It is also possible to derive estimators 
for non-zero total momenta, which is omitted here. 

Finally, one can derive an estimator for the isotope exponent of the polaron 
effective mass a^. In the (bi)polaron mechanism of superconductivity, is 
related to the isotope effect on the critical temperature j86] . The mass isotope 
exponent is defined as m* = M"^' , where M is the ion mass. Since the present 
method calculates the inverse polaron mass, is more conveniently expressed 
via the inverse effective mass 



M d ( I 



d 



fe) 



1 



(34) 



The last transformation follows from the scaling M cx The estimator for 
is derived directly from the estimator for the inverse effective mass A 
path's weight depends on the phonon frequency only via the polaron action 
g^poi^ (One should recall that the definition of a Monte Carlo average ((Z\r^)^) 
contains Apoi in the numerator and denominator.) Upon differentiation one 
obtains 



dA 



pol 



duj 



Ml, 



dA 



pol 



(35) 



In taking the frequency derivative, one should keep the combination Mui'^ 
constant, as it is independent of the ion mass. As with the phonon action 
and potential energy, the estimators for the number of phonons and isotope 
exponent can be split into a double sum over path's segments. The expressions 
for the respective summands can be found in |36j . 



4 Polaron properties 

The focus of this section is going to be the polaron properties that have 
been obtained with the continuous-time path-integral Quantum Monte Carlo 
method. For each topic, additional details of the algorithm will be provided, 
if necessary, on top of the general description presented above. 

4.1 Long-range interaction produces light polarons 

The Holstein model [67] describes polarons in systems with well-screened 
short-range electron-phonon interaction, such as molecular crystals. It is in 
the short list of "main" polaron models. Due to its simplicity, the model 
has been extremely popular with the numerical community. Almost any new 
numerical method was first tried on the Holstein model, for which a lot of 
accurate results are available. Also, for a long time the Holstein model was 
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Fig. 3. Properties of the Holstein polaron on d = 1,2,3,4,6 hypercubic lattices. 
The filled and open symbols show data for the phonon frequencies tvjjjt = 0.4 and 
1.0 (d = 1), 0.8 and 2.0 (d = 2), 1.2 and 3.0 (d = 3), 1.6 and 4.0 (d = 4), and 2.4 
and 6.0 (d = 6), respectively. The coupling constant is defined as A = Ep/{zt) = 
/ {2Mu>^ zt) , where z = 2d is the number of nearest neighbours. The effective mass 
is measured in units of mo = h^/{2ta^). Note that because the phonon frequency 
increases with d, the curves in different dimensions have similar shapes except for 
small differences around the A of polaron formation. 
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regarded as a "generic" raodel for all narrow-band systems where the e-ph 
and lattice effects were equally important. The locality of the interaction was 
considered simply a matter of convenience, and the Holstein model itself was 
believed to contain all the qualitative features of lattice polarons. This turns 
out not to be the case. The main reason lies in the very sharp, exponential, 
dependence on the model parameters of the main polaron characteristic, the 
effective mass. Some properties of the Holstein polaron are presented in Fig.[3l 
Any approximation or model simplification that have little effect on the po- 
laron energy can have dramatic consequences for the mass, leading to false 
qualitative conclusions. In particular, the Holstein model is not adequate for 
complex oxides of transition metals such as the cuprates or manganites. In 
those materials the bare bands are narrow, ~ 1 — 2 eV, so the lattice effects are 
important. At the same time, a low density of free carriers cannot fully screen 
the Coulomb interaction between electrons and distant anions and cations, 
leading to a long-range e-ph interaction. 

A model of this kind was considered a long time ago by Eagles [87] . (The 
corresponding polaron was called there a "nearly-small polaron".) More re- 
cently, a long-range model with linearly polarized phonons was put forward in 
relation to high-temperature superconductors in [32]. The model is depicted 
in Fig. SJa). An electron hops within a rigid chain (in the one-dimensional 




-2-1 

Fig. 4. (a) One-dimensional model with a long-range electron-phonon interaction. 
The electron moves along the bottom chain of sites shown by full circles (denoted 
by n). The vibrating ions shown by open circles oscillate along the z-axis. The 
interaction is characterized by the z-projection of the Coulombic force (|36ll . (b) 
Spatial profile of function (|30|) . that characterizes the retarded interaction of 
different parts of the polaron path. The parameters are k — 1 and d — a. Negative 
$ is shown in order to draw analogy with a potential well. This shape of # should 
be contrasted with that of the Holstein model ^hoI = n^S^^r. A smooth <!> cannot 
localize the path as well as a sharp one, which results in a smaller effective mass in 
the case of a long-range model. 
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case) or plane (in the two-dimensional case) and interacts with ions that are 
positioned above the electron chain (or plane) and vibrate perpendicular to 
the chain (plane). Such a polarization was chosen to represent the c-polarized 
phonon modes that feature prominently in the cuprates |88| . The force func- 
tion was taken to be the Coulomb force projected on the 2;-axis 



2(n-m)2 



/m(n) = na = , . (36) 

|n — m|-^ [(n - rim) + a J-^' 

Here a is the lattice constant, d is the distance between the vibrating ions 
and the electron plane, and rim is the vector of the electron Wannier orbital 
in the same unit cell as ion m. The interaction strength is characterized by 
the parameter n. Alternatively, one can use a dimensionless parameter 

m 

which is a ratio of the polaron energy in the atomic limit Ep (at t — 0) to 
the kinetic energy of a free electron zt (at k = 0). Here z is the number of 
nearest lattice sites. The confinement of the polaron path is determined by the 
function (p, (|30p . whose profile is shown in Fig. |4Kb). (It should be compared 
with the Holstein confining function K^6nn'-) Due to loose confinement, the 
polaron path diffused farther than in the Holstein case, which is equivalent to 
a much lighter mass. Such a polaron was named the "small Frohlich polaron" 
(SFP) in 




Fig. 5. Effective masses of small Holstein and Frohlich polarons in units of mo — 
l{2ta^). (a) d = 1. (b) d = 2. (c) The mass ratio of the small Holstein and 
Frohlich polarons for several model parameters. The ratio scales exponentially with 
the coupling, and could exceed 100. 



The masses of SFP and SHP calculated with PIQMC are compared in 
Fig. [5] SFP is slightly heavier at small A but much lighter at large A. The ra- 
tio of the two masses is a non-monotonic function. This observation was later 
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confirmed by exact diagonalization ,52J and variational [5H1[S2] methods. The 
most interesting property of SFP is exponential reduction of mass relative to 
SHP, see Fig. [SJc). This effect is independent of the dimensionality because 
it originates solely from the long-ranginess of e-ph interaction. The mass re- 
duction can be very large. For example, in two dimensions at Tioj/t = 0.5 
and A = 1.2, toshp ~ 220 while msFP is only w 9. In [36], an additional 
screening factor was added to the force function ([36|) . All polaron properties 
smoothly interpolated between those of SHP and SFP as the screening radius 
was changed from zero to infinity. 

The physical reason for the small mass of SFP is easy to understand. The 
mass is determined by the overlap of the ionic wave functions before and 
after the electron hop. In the case of short-range e-ph interaction, the lattice 
deformation must relax all the way back to the equilibrium state after the 
electron leaves the site. Likewise, the lattice deformation at the new electron's 
location must form anew from the equilibrium state. As a result, the overlap 
integral is exponentially small. In the case of a long-range interaction, the 
deformation at the new location is partially pre-existent. The ions must move 
by a smaller distance, which results in a much larger overlap of the wave 
functions. That also produces an exponentially large mass, but with a smaller 
exponent. (The fact that the mass is still exponential in coupling justifies 
the term "small" in the polaron's name.) This effect can also be understood 
within the Lang-Firsov (LF) approximation for the polaron mass £9, ,90 that 
becomes exact in the limit of infinitely fast ions, hu) ^ t. According to LF, 
the polaron mass is given by m* = toq e^'^^-pl^^^ ^ where is the polaron shift 
from (|37p. The dimensionless parameter 



depends on the shape of the c-ph interaction. For the force function at 
d = a, 7 = 0.387 for the linear chain, 0.334 for the square lattice and 0.320 for 
the triangular lattice. Note also that a long-range e-ph interaction smoothens 
the polaron crossover, cf. Fig. O which makes the single exponential form be 
more applicable for the entire coupling range [32j . 

These results demonstrate that the Holstein model is an extreme model 
that predicts the largest mass for the given polaron energy. As such, the 
Holstein model is not adequate for ionic systems with poor screening where 
the e-ph interaction is long range. 

4.2 Enhancement of anisotropy by electron-phonon interaction 

Two different behaviours of the lattice polaron can lead to an interesting 
situation when the same particle is Holstein-like and Frohlich-like along two 
different lattice directions. A concrete model that features such properties 
was put forward in |33j . The model is an extension of the force function (|36p 



7 = 1 



En./m(0)/m(a) 



(38) 
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to three dimensions, see Fig. IHlJa). Clearly, the e-ph interaction is anisotropic, 
which translates to an even higher anisotropy of polaron spectrum. Upon hops 
along the z direction, the electron must reverse the sign of deformation. Thus 
the ion overlap integrals are even smaller than in the Holstein case. In contrast, 
the movement in the xy plane is Frohlich-like, as described in the previous 
section. The deformation is partially prepared, the overlap integrals are large, 
and the effective mass is small. This reasoning is corroborated by the spatial 
profiles of the path confinement function ^, see Fig. IHJb). When the path dif- 
fuses along the x direction, ^ is smooth which translates into easy "escape" . 
Along the z direction, is much steeper. In fact it even changes sign upon 
migration to the nearest plane, which reflects the sign reversal of the lattice 
deformation after the electron hop. Since affects the path's weight exponen- 
tially, one expects exponentially large difference between the polaron masses 
in z and x direction: m*/m* oc e^'^--'^-)^/'^, where A = I.^iZk^ I{Y1MujH^) 
and Lo = hiojix is the dimensionless phonon frequency. As a result, the mass 
anisotropy grows exponentially with A and can reach very big values. In addi- 
tion, the mass anisotropy is a sharp function of the phonon frequency, which 
should lead, for example, to an isotope effect on anisotropy. 
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Fig. 6. (a) A three-dimensional polaron model with anisotropic e-ph interaction. 
The electron moves within and between the planes of filled circles with hopping 
integrals —t^ and — tz, respectively. It interacts with vibrating ions shown by open 
circles. The vibrations are polarized along the z-axis. In the case of a long-range 
interaction the polaron is Frohlich-like in the xy plane and Holstein-like along the 
z direction, (b) The (negative) polaron path confinement function —K~^<P{x,y,z). 
An additional factor e"'"""*'/^ with R = 10 lattice constants in the x direction, 
has been added to the force function (|36p to help the lattice sum to converge. The 
confining profile along z is much steeper than along x. 
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The results of PIQMC analysis of the model are summarized in Fig. [7] 
[55] , The inset shows a typical behaviour of the polaron masses at Tiuj — 1.0 1^, 
and tz = 0.25 tx- (This choice of hoppings ensures the isotropy of the bare 
spectrum, since the lattice constant in z direction is assumed to be twice 
the one in x direction.) As expected, m* grows exponentially with coupling. 
Interestingly, the z mass grows sM|jer-exponentially with a large quadratic 
component: m* « nixo e^-'^^^^^-^^'^ . It is very possible though that this is 
still transient regime, and to* approaches pure exponential growth at still 
larger A (at which becomes so large it is difficult to calculate). The mass 
anisotropy for several sets of model parameters is shown in the main panel 
of Fig. [71 Due to the super-exponential growth of m*, the anisotropy is also 
super-exponential, e.g. TO*/m*y « g0.42A-i-o.7iA ^^j. _ 1.0 (circles). At a 
smaller frequency hu = 0.5 (squares) the anisotropy is exponentially larger, 
as expected from the reasoning given above. This implies the existence of an 
isotope effect on the mass anisotropy. The third model parameter, the bare 
hopping anisotropy, turns out to have little effect on the anisotropy of the 
polaron spectrum. The mass anisotropy for tz = 0.1 1^ and the same phonon 
frequency Uui — 0.5 tx is shown by triangles in Fig. [T] While being 2.0 — 2.5 
times larger at small A, the anisotropy approaches that of the t^ = 0.25 tx 
case at large A. It shows that it is primarily the e-ph interaction that governs 
the polaron anisotropy. This conclusion is further supported by studies of the 
Holstein model with anisotropic bare hopping [33j where no enhancement of 
the polaron anisotropy was observed. 
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In relation to the effect described, it is interesting to note a well-docu- 
mented discrepancy between the theoretical and experimental anisotropy of 
the cuprates [HI] ■ According to band structure calculations, the anisotropy of 
resistivity of LSCO and YBCO compounds should be about 10-30 [5^155] . 
At the same time, the experimental anisotropy of resistivity is between 10^ 
and 10^ depending on the level of doping. The anisotropy of bismuthates is 
even higher, (5 — 80) ■ 10''' [91], which is difficult to explain with the conven- 
tional Bloch-Boltzmann theory. According to the present results, anisotropic 
interaction with z-polarized phonons is a sufficient condition for a very large z 
effective mass. Of course, the anisotropy of mass and resistivity are two differ- 
ent things, but it would be fair to assume that the former is at least partially 
responsible for the latter (see, e.g., [91] page 85). This idea still awaits proper 
development. Alternative explanations of the anomalous z-transport exist as 
well [M]. 

4.3 Spectrum flattening and polaron density of states 

A great power of the PIQMC method is the ability to calculate an entire po- 
laron spectrum in any dimensionality in a single run, as was explained in Sec- 
tion [331 This opens up an exciting possibility to calculate the polaron density 
of states (DOS) p[E) = N''^ X^k '^(■^ - ^^k + £^0) by discretizing the energy 
interval and histogramming E^i values at the end of simulations. The polaron 
spectrum in the adiabatic limit (t ^ Tllo) possesses an interesting property of 
flattening at large lattice momenta, which has been well documented in the 
literature [50, 5l]|55l[95]. In the weak-coupling limit, the flattening can be 
readily understood as hybridization between the bare electron spectrum and 
a momentum-independent phonon mode. The resulting polaron dispersion is 
cosine-like at small K and flat at large K. As larger couplings this picture 
becomes less and less intuitive but the flattening remains as a matter of fact 
|96j . It is noteworthy that in large spatial dimensions the phase space is dom- 
inated by states with large momenta. If all those states have close energies, 
they should form a peak in the DOS close to the top of the polaron band. Ex- 
act PIQMC calculations have confirmed that this is indeed the case f30]. The 
evolution of the DOS of the isotropic Holstcin model with phonon frequency 
u) in two and three dimensions is shown in Fig.[8l^a) and (b). In all the cases 
presented the polaron is fully developed, with the bandwidth much smaller 
than u>. 

At small DOS develops a massive peak at the top of the band. The 
peak is more pronounced in d = 3 than in d = 2. The van Hove singularities 
are absorbed in the peak and as such cannot be seen. With increasing w, the 
polaron spectrum approached the cosine-like shape in full accordance with 
the Lang-Firsov non-adiabatic formula. The respective DOS gradually assume 
the familiar shapes of tight-binding zones with renormalizcd hopping integrals. 
The van Hove singularities are clearly visible. These results have an interesting 
corollary for the Holstein model. At small-to-moderate lo in two and three 
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Fig. 8. (a) and (b) The evolution of the density of states of the Holstein polaron 
with phonon frequency a; in d = 2 and 3, respectively. Each graph was obtained by 
calculating the polaron spectrum at 100,000 K points randomly distributed in the 
Brillouin zone, and histogramming the results between 50 energy bins. Each spec- 
trum point was calculated by averaging 250,000 values of cosKZir taken every 10th 
path update. Every 5000 measurements the path was reset and then equilibrated 
for 1000 updates, (c) The same for the small Frohlich polaron in d = 2. 



dimensions, the bottom half of the polaron band contains a tiny minority of 
the total number of states (measured in low percentage points). In most real 
systems those states will be localized and hence irrelevant. All the system's 
responses will be dominated by the states in the peak. The physical properties 
of these physically relevant states are going to be quite different from those 
of the ground state that are usually analyzed theoretically. 

It is interesting to look at the DOS of the long-range model ([36]). The two 
dimensional DOS is shown in Fig.[51^c). It is much closer to the tight-binding 
shape than the Holstein DOS at the same parameters. The polaron spectrum 
and density of states is another manifestation of the extremity of the Holstein 
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model. Long-range e-ph interactions remove those peculiarities and make the 
polaron bands more "normal" . 

Some densities of states of the Holstein model with anisotropic hopping 
were presented in [30j . 

4.4 Isotope exponents 

Isotope substitution is a powerful tool in determining if a particular phe- 
nomenon or feature is of phonon origin. Some polaron properties, such as 
mass, are sharp functions of the lattice parameters, therefore isotope effects 
in (bi)polarons are strongly enhanced. In particular, anything that depends 
on the (bi)polaron mass should exhibit a large isotope exponent [53] . A large 
isotope effect was observed, for example, on the magnetic penetration depth 
in cuprates [HZl [551 [SS]j which was interpreted as evidence for bipolaronic 
carriers in the superconducting state. 

As was explained in Section [?751 the PIQMC method enables approxima- 
tion-free calculation of the isotope exponent on the (bi) polaron effective mass 
for a large class of e-ph models. A good feel for a in the strong-coupling 
limit can be obtained from the anti-adiabatic expression for the polaron 
mass, TO* = Too e''^-^p^^^K Since the polaron shift Ep is ion-mass-independent, 




Fig. 9. Mass isotope exponents of the small Holstein and small Frohlich polarons 
in different dimensions and for different phonon frequencies. The thin solid lines 
indicate the strong-coupling limit a = ■ Note that a of the Holstein polaron in 
d = 2 and 3 is negative at small cu. 
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" ~ 2 huj' Thus a is proportional to the couphng constant A. In the weak 
coupUng regime, a can be computed perturbatively. For example, for the one- 
dimensional Holstein model, the second order yields 

^(2) ^ u^u^ + 2u + A) 

%dHolstcin = ^ , ,-.5 , (39) 

where uj = hcu/t. (The second-order coefficients of other polaron properties of 
the one-dimensional Holstein model can be found, e.g., in |36| . The second or- 
der coefficients for higher-dimensional lattices were given in '37] . ) The isotope 
exponent is again proportional to A but with a different coefficient. The two 
linear dependencies should be smoothly connected over the polaron crossover. 

The mass isotope exponents for the small Holstein polaron are shown in 
the top row of Fig. [SI Before the polaron transition a is small reflecting the 
non-polaronic character of the carrier. Note that in d = 2 and 3, a is negative 
at small frequencies After the polaron crossover, which always begins 
at A ~ 1 but ends at A that increases with uj, the isotope exponent quickly 
reaches the strong-coupling asymptotic behaviour with 7 = 1. Notice that 
the beginning and the end of the transition are clearly identifiable on the 
plots. Thus the mass isotope exponent can be useful in analyzing the Holstein 
polaron transition. 

The case of the small Frohlich polaron is somewhat different, see the bot- 
tom row of Fig. [9l At small coupling, the exponent grows linearly with a 7 
close to the strong coupling limit, but then deviates to smaller values. The 
exponent returns back to the strong-coupling limit at much larger A, whose 
value decreases with increasing uj. This final approach happens when the en- 
tire path is mostly confined to one lattice site, and only rarely deviates to the 
first nearest neighbour. This interesting behaviour is not yet fully understood. 

As was shown in the previous section, in the adiabatic regime the prop- 
erties of the polaron mass do not fully represent those of the vast majority 
of polaron states. Additional insight can be gained from an isotope effect on 
the polaron bandwidth or even on individual spectrum points [35j. The iso- 
tope effect on polaron spectrum and density of states in d = 2 is illustrated 
in Fig[TOl The ratio of the two phonon frequencies, fiu; — 0.80 i and 0.75 
has been chosen to roughly correspond to the substitution of ^^O for in 
complex oxides. One can see that the polaron band shrinks significantly, by 
20-30%, for both polaron types. The middle panels show the isotope exponents 
on spectrum points calculated as 

1 (uj) AEk , , 

where the angular brackets denote the mean value of either the two frequencies 
or of the two energy values. An interesting observation is that aK of the 
Frohlich polaron is roughly independent of K (±10%). In the Holstein case 
Qk dips in the vicinity of the F point. (Small fluctuations at intermediate 
momenta are due to statistical errors in averaging (cosKZir).) 
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Fig. 10. Isotope effect on polaron spectrum and density of states. Top row: the 2d 
Holstein polaron at A = 1.2. Left: polaron spectrum at tj = 0.80 and 0.75. Middle: 
the isotope exponent for each K-point. Right: the density of states for the two 
frequencies. Bottom row: the same for the small 2d Frohlich polaron at A = 2.4. 

4.5 Jahn- Teller polaron 

Density-displacement is only one possible type of e-ph interaction. Examples 
of other types are the deformation potential and Su-Schrieffer-Heeger (SSH) 
interaction. In the former, the electron density interacts with a gradient of 
lattice displacement whereas in the SSH case the lattice deformation is cou- 
pled to electron's kinetic energy. The Jahn- Teller ( JT) interaction is one of the 
most complex types because it usually involves a multidimensional electron 
basis and a multidimensional representation of the deformation group. The JT 
interaction is active in some molecules and crystals of high point symmetry. 
The JT effect was also a guiding principle in the search for high-temperature 
superconductivity [lOOj . More recently, a JT theory of the cuprates was devel- 
oped in [101[|102lll03| . There exist different flavors of the JT interaction. The 
simplest one is the E ® e interaction [104j that describes a short-range cou- 
pling between twice-degenerate Cg electronic levels (ci, C2) and a local double- 
degenerate vibron mode {Ct])- The Hamiltonian reads 

(nn') 

- X! (clilCn2 + c]^2Cnl) ?/n + (cj^lCnl - c|,2Cn2j Cn 
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2M \dQ ^ drjl 



(41) 



The symmetry of the interaction ensures the same couphng parameter k for 
the two phonons. The kinetic energy is chosen to connect Uke electron orbitals 
of the nearest neighbours. This choice is somewhat arbitrary (makes it more 
"Holstein-like" ) and not dictated by symmetry. Because the ionic coordinates 
of different cells are not coupled, the model describes a collection of separate 
clusters that are linked only by electron hopping. To relate the Hamiltonian 
to more realistic situations, phonon dispersion must be added [1011 1105] . 

An important property of the E ® e interaction is the absence of an exact 
analytical solution in the atomic limit t = 0. In the density-displacement case 
the dynamics of each C, is described by a one-dimensional differential equation 
under a constant force. That has a shifted oscillator solution, which serves 
as a convenient starting point for various strong-coupling expansions. Here, 
in contrast, the atomic limit is described by two coupled partial differential 
equations for the electron doublet Vi,2(Cj^)- Although it is possible to sepa- 
rate variables and reduce the system to two ordinary second-order equations, 
they seem to be too complex to admit an analytical solution. At large cou- 
plings, however, the elastic energy assumes the Mexican hat shape and the 
phonon dynamics separates into radial oscillatory motion and azimythal ro- 
tary motion. This results in an additional pre-exponential factor oc k in the 



ion overlap integral, leading to the effective mass m*^^^ = mo 



where 



9' = j^m- 

A path integral approach to Hamiltonian (|4T|) was developed in f34]. Be- 
cause there are two electron orbitals, the electron path must be assigned an 
additional orbital index (or colour) a — 1,2. Colour 1 (or 2) of a given path 
segment means that it resides in the first (second) atomic orbital of the elec- 
tron doublet. The short-term density matrix is 



p{AT) = {r'a';{aWn}\e 



-AtH 



ra; {Cn}{??n}) 



,(42) 



v4ph(Z\T) = {KAT){Sal ~ (5a2)Cn 

- E {i^2^[(Cn - C;)^ + (^n - V'J'] + iAr)^{Cl + vl)]- (43) 

Here a — 1(2) when a — 2(1), that is a is "not" a. This expression reveals a 
difference between the two phonons. Phonon C is coupled to electron density, 
like in the Holstein case. The difference from the Holstein is that the direction 
of the force changes to the opposite when the electron changes orbitals. In 
contrast, phonon ry is coupled to orbital changes themselves: the more often 
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the electron changes orbitals, the more "active" is rj. (Discrete orbital changes 
are analogous to electron hops between discrete lattice sites, and as such are 
associated with "kinetic orbital energy" . Phonon coupling to orbital changes 
is analogous to phonon coupling to electron hopping in the Su-SchriefFer- 
Heeger interaction |107| .) Multiplication of a large number of the short-time 
propagators generates multiple terms, each of which contains a finite number 
of lattice hops and orbital changes. On top of that, the exponents combine in a 
total phonon action Ap\i{P) that comprises the free r] action Ar^o and the action 
of C phonons under an external force [A(; is similar to ipT|) ]. The next step 
is integration over the paths r]{T) and C(''")- Integration over C(t) is standard 
and performed as described as for the density-displacement interaction. The 
result is a factor e^'^ in the path's statistical weight, where 

A^[r{T),a{T)] =K^ / drdT'G(T - r') [6a(r).a{T'} - Sa{r)Mr')] ' (^4) 

Jo Jo 



G{r-r') 



2Muj 



-hujW- 



5r(T),r(T') + e ^'^ ' dr(r),r(T')+sgn(T-T')Zir 

(45) 

The last expression is valid under the condition e^^^"^ 1. As expected, action 
Aq is an explicit functional of the spatial path v{t) and of the orbital path a(r). 
Note that the C, phonon favours like orbitals and disfavours orbital changes. 
Integration over "qij) is trickier for it contains 77 as pre-exponential factors. A 
typical term with Ns orbital changes has the following form: 

(«Zir)"ryr(.,)(ri)r7r(r.)(r2)...r;r(.„j(T^Je^''«[""Ml . (46) 

Here, Ts is the time of sth orbital change (s — 1, . . . , A^^), y{ts) is the electron 
position at this time, and ?7i.(ts)(t's) is the 7y-displacement at this site at this 
time. For odd A^^, path integration over rjnij) produces zero by parity. For 
even m, the integration can be performed by introducing fictitious sources, 
calculating the generating functional and differentiating it m times. The result 
is a sum of all possible pairings of Tj. Within each term, each pair contributes 
the factor k^G{ts—Ts' ) , with G given by . Since G is positive, the statistical 
weight increases with the number of orbital flips. Thus, the 77-phonon favours 
orbital flips, in complete contrast with the C-phonon. Since both phonons are 
governed by the same Green's function (which should be no surprise since 
the two phonons are related by symmetry) one should expect a dynamical 
equilibrium between the two phonons and some finite mean value of orbital 
flips for given coupling. The shift partition function is 



^4r = V E E / ••• / (dT)^Mdr)^=W^M, (47) 

Nk=0,l,... N,=0,2,...^ •'^ 



{ pa 



g-4c[r(r),a(r)] ^ ^^g-j 
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Compared to (|27p . the above expression involves additional multiple integra- 
tion over the times of orbital flips. Each term is associated with a spatial- 
orbital path with Nk kinks and (Ng/2) ?7-phonon pairings, as shown in 
Fig. lllf a). A Markov process is organized by inserting/removing spatial kinks 
and by attaching/removing the pairing lines. The acceptance rules for the pair- 
ing lines can be derived by extending the method of Section 13.11 to two-time 
updates. For details see [34]. The path shown in Fig. [TlT a) can be considered 
a space-time diagram. In fact, the update process just described is a version 
of the Diagrammatic Monte Carlo method [53] • 




Fig. 11. Extensions of the basic polaron PIQMC method, (a) In the Jahn- Teller 
polaron, the path carries an internal orbital index shown as black or white sections of 
the path. The r;-phonon parings are represented by dashed lines. Each vertex changes 
the index, (b) In the bipolaron case, the system is represented by two imaginary-time 
paths. A direct configuration is shown on the left and an exchange one on the right. 
The interaction of segments marked A or B contribute to the energy of individual 
polarons. The segments that belong to different paths, such as those marked C, 
contribute to the interaction between the polarons. 



After the update rules are established, the JT polaron properties can be 
calculated with no approximations. The mass, spectrum, and density of states 
are obtained as for the conventional lattice polarons. For the JT energy the 
thermodynamic estimator yields: 

\ ' shift 

The first two terms are familiar from the Holstein case. The first is the kinetic 
energy of the electron. The second is the potential energy of interaction with 
the C plus the (positive) elastic energy of C^. The other two terms are new. 
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The third one represents the negative "orbital energy" caused by orbital flips. 
These processes excite phonons rj, which results in the positive fourth term. 

The number of excited phonons is an important characteristic because 
it provides an internal consistency check for the algorithm: one expects equal 
mean number of C and 77 phonons. The phonon number estimator is derived as 
the derivative of the free energy with respect to {Hlo) under a fixed combination 
/{Muj). One obtains 



dAc 



ph,C 



d{huj) 



(50) 



shift 




' shift 

Results of PIQMC calculations are shown in Fig. [12] [S^- Most properties 
behave similarly to those of the d = 3 Holstein polaron at the same phonon 
frequency, huj = 1.0 i. For example, the kinetic energy. Fig. [TWa). sharply 
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Fig. 12. Physical properties of the d = 3 Jahn- Teller polaron at Tuu — 1.0 1 [34) . 
(a) The total and kinetic energy, (b) The effective mass compared with the d = 3 
Holstein polaron at the same phonon frequency, (c) The number of excited phonons 
of both types, (d) The density of states of the JT polaron at A = 1.3. 
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decreases by absolute value between A ~ 1.2 and 1.4. The JT polaron mass is 
slightly larger at the small to intermediate coupling, but several times smaller 
at the strong coupling. This non-monotonic behaviour of the ratio of the JT 
and Holstein masses was later confirmed by accurate variational calculations 
[60] . although in that work the JT polaron (and bipolaron) was investigated 
in one spatial dimension. The relative lightness of the JT polaron is consistent 
with Takada's result mentioned above , 106j . The number of excited phonons 
of both types is shown in Fig. [T2tc). As expected, numerical values coincide 
within the statistical error, which validates the numerical algorithm. Interest- 
ingly, the shape of the phonon curves is similar to that of the logarithm of the 
effective mass. This suggests an intimate relationship between the two quan- 
tities, again similarly to the Holstein case. Finally, the density of JT polaron 
states features the same peak at the top of the band, caused by the spectrum 
flattening at large polaron momenta. 

In summary, the locality of the JT interaction and the independence of 
vibrating clusters result in the same extremity of polaron properties as in 
the 3ci Holstein model. One should expect that either a long-range JT or 
phonon dispersion will soften the sharp polaron features and make them more 
"normal" . This is an interesting research opportunity for the PIQMC method. 

5 Prospects 

In this Chapter, a theoretical approach to lattice polarons based on the statis- 
tical path integral has been developed. A characteristic feature of the method 
is a tight integration of analytical and numerical elements within one tech- 
nique. The numerical part of the algorithm is greatly aided by the preliminary 
analytical steps: path integration over the ionic coordinates and reduction of 
the double integral over imaginary time to a discrete sum over path segments. 
Thus, on the way from the Hamiltonian to observables, exact analytical trans- 
formations are carried as far as possible increasing the accuracy and speed of 
the simulations. 

The PIQMC polaron method is still very much under development. Due to 
the space constraints, only key elements have been presented here and many 
technical details have been left out. This concluding section will be used to 
outline various extensions of the method. Some of them have already been 
done, others can be done but await proper implementation. Certain problems 
have not yet been solved even on the conceptual level, and those will be 
discussed too. 

5.1 Bipolaron 

The possibility of two polarons forming bound pairs, or bipolarons, opens a 
whole new dimension in polaron physics |108[ 1109] . The two major competing 
factors are the phonon-mediated attraction between the polarons that favours 
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paring and direct Coulombic repulsion that prevents pairing. Kinetic energy, 
exchange effects and short-range lattice effects play important supporting roles 
that can tip the balance toward or against paring. All of this may lead to a rich 
and complex phase diagram. The bipolaron concept has been well-developed in 
the literature. The large continuous-space (Frohlich) bipolaron in ionic media 
was considered in [201 1110| llllj - and the small, localized, lattice bipolaron 
in |112j . Since the bipolarons are charged bosons at low density n, they can 
superconduct below the Bose-Einstein condensation temperature 

Tbe = 3.31- -. 52 

Kb ni** 

The concept of mobile lattice bipolarons and bipolaronic superconductivity 
was put forward by Alexandrov and Ranninger in 1981 [113) and since then 
thoroughly developed by Alexandrov and co-workers [701 Ell IHSl [Ml 11141 11151 
111611117] . Note that Tbe is inversely proportional to the bipolaron mass. This 
suggests the following recipe to increase the critical temperature according to 
the bipolaronic mechanism: reduce the effective mass while keeping the bipo- 
larons non-overlapping. (Clearly, simply reducing the mass will increase the 
bipolaron radius and reduce the overlap density. There is a trade-off between 
the n^/"^ and (m**)^^ factors in ([5^ . which creates a challenging optimization 
problem.) 

As with polarons, the first application of path integrals to the bipolaron 
was to the continuum-space Frohlich bipolaron IllOj . Until very recently, 
the only application of path integrals to the lattice bipolaron was the work 
by de Raedt and Lagendijk on the discrete-time version of their path- 
integral QMC method. The authors limited the consideration to the extreme 
adiabatic limit and computed the phase diagram of the Holstein bipolaron 
with an additional Hubbard repulsion. A continuous-time version of PIQMC 
was developed very recently and the first results presented in [3^ . 

The lattice bipolaron, mostly with Holstein interaction, has been studied 
by a variety of other numerical methods: exact diagonalization [5Dllll8j . vari- 
ational [57l [58l [60] , density matrix renormalization group [66] , diagrammatic 
quantum Monte Carlo [H] , and Lang-Firsov quantum Monte Carlo [13 1119j . 

Conceptually, generalizing the PIQMC method to the bipolaron is straight- 
forward. There are two paths that visualize the imaginary-time evolution of 
the two fermions. To enable calculation of the effective mass, open bound- 
ary conditions in imaginary time must be used: the top ends of the paths 
must be the same as the bottom ones up to an arbitrary lattice translation 
Ar. Singlet and triplet states can be separated by allowing the paths to ex- 
change, see Fig. lll( b). Then the bipolaron spectrum, mass, and singlet-triplet 
splitting is calculated as explained in Section [21 Phonon integration is per- 
formed as before with the same resulting action, ([^^ -(|26 p . The difference is 
the force that is acting on the oscillator m at time r. In a bipolaron system 
it receives contribution from both electrons, fmir) — /m[i"i(T)] -|- /m[r2(T)]. 
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Since the action is quadratic in /, it becomes a sum of three terms that con- 



tain the products /[ri(r)]/[ri(r')], /[r2(r)]/[r2(T')], and /[ri(T)]/[r2(r')], 



respectively. The first two terms describe the retarded self-interaction of the 
two paths and are responsible for polaron formation. The third cross-term 
describes the interaction between the paths. It is responsible for the interac- 
tion between the polarons and for bipolaron formation. This is illustrated in 
Fig. IllT b). The estimators for various observables and specific expressions for 
the segment-to-segment contributions are derived analogously to the polaron 
case. Technical details of the derivation have not yet been published, and will 
be published elsewhere. 

One interesting effect that has already been considered by the PIQMC 
method is the "crab-like" bipolaron that can exist on certain lattices [39j . 
Usually, the bipolaron is much heavier that the constituent polarons, its mass 
scaling as the second power of the polaron mass in the non-adiabatic regime 
and as the fourth power in the adiabatic limit [120j . (Path integrals provide a 
useful visualization of this effect, see Fig. [TlT b). At a small phonon frequency, 
the two paths interact over their entire length. Therefore they are much more 
difficult to separate, which results in slower imaginary-time diffusion, and 
hence a heavier particle.) However, on the triangular, face-centered cubic and 
some other lattices the intersite bipolaron can move without breaking. This 
results in an effective mass that scales only linearly with the polaron mass. 
This effect was predicted some time ago by Alexandrov [70], and recently 
confirmed by exact PIQMC simulations in [ 39] . 




5.2 Further extentions 

Recall that the basic Hamiltonian (fT^ contains the simplest possible form 
of the electron kinetic energy and the simplest possible model of the crystal 
lattice. It will now be explained how to relax these restrictions. First of all, 
it is straightforward to add electron hopping beyond the nearest neighbours. 
Because the path is defined in real space, this only requires the introduc- 
tion of additional kink types that move the path in the respective directions. 
The values of the distant hopping integrals can be absolutely independent of 
the nearest-neighbour values as long as they are negative. The methods of 
Section 13.11 is readily generalized to paths containing kinks of different mag- 
nitude. The statistical weight p5]) is now a product of factors (tiAr) along 
the path. The stochastic acceptance rules now contain the factor {tiP) where 
ti is the hopping amplitude of the proposed kink. It is interesting though that 
contribution to the kinetic energy is the same (equal to for all kinks inde- 
pendently of the ti value. It is the mean number of kinks that is affected by ti, 
leading to different partial kinetic energies along different lattice directions. In 
exactly the same manner one can study models with anisotropic (but nearest- 
neighbour) electron hopping. This was done in [301 133], cf- Section There 
remains one important restriction on the values of the hopping amplitudes: 
they all must be negative to ensure the positivity of the path's statistical 
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weight. If even a single ti is positive, some paths will continue being posi- 
tive but some will have negative weights. That will introduce a sign-problem, 
which will eventually render the simulation numerically unstable. Another ex- 
tension related to the kinetic energy is lattices of different symmetries. Again, 
all that is changed is the table of kink types that specifies which two spa- 
tial points (lattice sites) each particular kinks connects. In this way, PIQMC 
has been applied to the triangular, face-centred-cubic, and hexagonal Bravais 
lattices in [37] . 

It should be mentioned that PIQMC can also simulate the (bi)polaron in 
dimensionalities larger than 3, see Fig.[3l (The d = 4 Holstein polaron was pre- 
viously investigated in [59].) Although such calculations have little practical 
value, they can still be useful in assessing the accuracy of approximate meth- 
ods that rely on the number of spatial dimensions being large. (An example 
of such a method would be the Dynamical Mean Field Theory [121] .) 

Consider now the phonon subsystem. It is well-known that a bosonic path 
integral can be calculated analytically for any quadratic action. This means 
that the ionic coordinates can be eliminated even when they are coupled, 
that is in the case of dispersive phonons. Moreover, integration can be done 
for any number of degrees of freedom per unit cell. De Raedt and Lagendijk 
studied a one-dimensional Holstein polaron with phonon dispersion as early 
as in 1984 [2j. They found that the critical coupling of polaron formation 
decreases as the phonon mode becomes soft. To our knowledge, this remains 
the only exact analysis of a model with dispersive phonons to-date. A general 
phonon integration for the purposes of PIQMC method was performed in [38] . 
Conceptually, the result is similar to P^ - ([?7|l . The action is a sum of a pe- 
riodic contribution and a correction due to the shifted boundary conditions. 
The action is a double integral over the imaginary time. However, there is an 
additional sum over the phonon momentum and branches. In addition, the ac- 
tion involves the eigenvalues and eigenvectors of the dynamical matrix. Thus, 
the polaron action comprises full information about the phonon spectrum and 
polarizations. 

An important advantage of the PIQMC method is the ability to study 
temperature effects. At a finite temperature T, the projected partition func- 
tion Zk receives contributions from high-energy states with momentum K. It 
is therefore not meaningful to calculate the polaron mass and spectrum. In- 
stead, the conventional periodic boundary conditions in imaginary time must 
be used, which makes all states contribute to the full thermodynamic partition 
function Z. The parameter Ar is set to zero. The remaining periodic action 
(|23p is valid at any temperature, large or small. This enables approximation 
free calculation of thermodynamic properties of the polaron: the internal en- 
ergy, specific heat, number of excited phonons, and static e-ph correlation 
functions. 

The grand challenge for PIQMC, as for almost any Quantum Monte Carlo 
method is efficient simulation of many-body systems. If the number of particles 
is three or more, the fermionic sign-problem is in general unavoidable. It is 
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possible that the sign-problem in many-ferniion e-ph models will be less severe 
that in the case of repulsive electron-electron models. Consider the limit of 
low electron density and strong e-ph interaction. The polarons will be bound 
in bipolarons, and the paths will fluctuate in pairs. To account for statistics, 
the top path ends must be constantly permuted, with even/odd permutations 
contributing -1-1/ — 1 phase factor to the partition function. However, any odd 
permutation will have to separate a path pair. This will increase the energy 
of the configuration and such an update will likely be rejected. In contrast, 
some even permutations will exchange whole bipolarons, which will keep the 
paths together with no significant energy increase. Such updates will likely be 
accepted. As a result, there will be many more even than odd permutations, 
and the average sign will be well defined. The stronger the coupling the more 
pronounced this effect will be, and the more the system will approach a purely 
bosonic limit (which is sign-problem- free). 

The temperature-control capability offers another unique opportunity: cal- 
culation of the Meissner fraction and superconducting critical temperature. 
The method was devised by Scalapino et al |122j and is based on the sum 
rule for the off-diagonal current-current correlation function. By comparing 
its static limit with the kinetic energy, one could determine a temperature 
at which the two quantities start to differ. This signifies the appearance of a 
Meissner fraction and determines the Tc. 

Finally, we comment on the calculation of dynamical (bi)polaron prop- 
erties such as the spectral function or optical conductivity. This is a classic 
and difficult problem for most existing Monte Carlo methods. Derivation of 
real-time correlators from imaginary-time correlators amounts to solving an 
ill-posed inversion problem, which is usually regularized by methods such as 
maximum entropy [123] . Fade approximants, or stochastic optimization [8]. 
Any of this techniques can be used in conjunction with the FIQMC method. 

5.3 Conclusions 

In summary, path integrals have played and continue to play an important 
role in the development of polaron physics. The path- integral quantum Monte 
Carlo method has proven to be a powerful and versatile tool. It has some 
unique advantages, such as system-size independence, the ability to calculate 
the density of states, mass isotope exponents, and temperature effects. Com- 
pared to the exact diagonalization, density-matrix renormalization group, or 
variational methods, FIQMC has a larger statistical error, but provides unbi- 
ased estimates for the (bi)polaron properties. Several novel qualitative results 
have been obtained with the FIQMC method: small polaron mass in models 
with long-range electron-phonon interaction, enhancement of the anisotropy 
of the polaron spectrum, a peak in the polaron density of states in the adia- 
batic limit, and the isotope effect on the polaron spectrum. These and some 
other effects have been discussed in this Chapter in detail. 
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